We use three different shapefiles for the continental U.S. land mass, the State waters of maine, new hampshire, massachusets, connecticut, rhode island, new york, new jersey, delaware, maryland, virginia, north carolina and the North East EEZ.
1.- Land shapefile; covers the US land territory for visualization. Data provided from the map package.
2.- State waters; covers the state waters of the NE US states. Data from data.gov.
License: No license information was provided. If this work was prepared by an officer or employee of the United States government as part of that person’s official duties it is considered a U.S. Government Work.
3.- EEZ shapefile; Used the Sea Around Us shapefle updated to June 2016.
unique(land_sf$ID)
As a first step we need to divide the NE US EEZ among the different states. For that we expanded state waters up to the 200 nautical miles to then estimate the percentage that each expanded-state-waters occupied. Note that in all cases these areas overlapped and percentages accounted for it. We did this by following these steps:
We set a buffer of 410000 m (410 km, ~ 221 nm) that overshoots the EEZ a bit, but is eventually cropped
# Buffer state waters
state_bf = st_buffer(state_sf, 410000) %>%
st_transform(4326) %>% # to match shape
st_set_crs(4326)
# ------------------------------ #
# Testing and visualizing the buffer
# ------------------------------ #
state_bf_plot <- ggplot() +
geom_sf(data = eez_sf, aes(), fill = NA) +
geom_sf(data = state_bf, aes(color = state), fill = NA)
#
# ggsave("../Results/Partial/state_waters_buffer.png",state_bf_plot)
# ------------------------------ #
state_bf_plot
Here we expand a grid within the buffer so we can estimate the proportion of each state
Once we have a gridded area, we converted the grid to a sf so we can merge it with the buffered states and finally filter out everything outside the states polygon
# ---------------- #
# Convert grid to sf
# ---------------- #
grid_sf <- st_as_sf(ne_grid,
coords = c("lon", "lat"),
crs = 4326) %>%
st_join(state_bf) %>%
filter(!is.na(state))
# Create data frame for future computations
# Note, will be used in next chunk
grid_bf_dt <- as.data.frame(grid_sf) %>%
select(index,state)
# group_by(state) %>%
# summarise(length(index))
# ---------------- #
# [TEST]
# Visualization of grid
# ---------------- #
gridExtra::grid.arrange(
# Overall (overlapping) position
ggplot(grid_sf) +
geom_sf(aes(color = state), alpha = 0.3) +
geom_sf(data = eez_sf, aes(),fill =NA) +
theme(legend.position = ""),
# Showing each state separatley
ggplot() +
geom_sf(data = grid_sf, aes(color = state),size = 0.1, alpha = 0.3) +
facet_wrap(~state) +
theme(legend.position = "top"),
nrow = 1)
Finally, we crop the grided buffers to within the EEZ to capture the actual water space.
Note: This step takes quite a while because of the size of the EEZ shapefile. No, you cannot use st_simplify() here
Here we interpolate the survey data within the previously created grid.
wtcpue < 0This is the main function used to interpolate the data per year. It follows a Triangular Irregular Surface method using the interp::interp() function.
tis <- function(input_data, grid, yr, taxa, reg){
# --------------- #
# Testing
# print(paste(yr))
# yr = 1976
# --------------- #
# Filter data
data <- input_data %>%
filter(year == yr,
spp == taxa,
region == reg
) %>%
# Average duplicated hauls in the same spot
group_by(region,year,spp,lat,lon) %>%
summarise_at(vars(wtcpue),
mean,na.rm = T) %>%
filter(wtcpue > 0)
# Only interpolate cases where there is more than 3 rows
# Marked by the function
if(nrow(data) <= 3){
fit_tin <- tibble()
}else{
# Triangular Irregular Surface
fit_tin <- interp::interp( # using {interp}
x = data$lon, # the function actually accepts coordinate vectors
y = data$lat,
z = data$wtcpue,
xo = grid$lon, # here we already define the target grid
yo = grid$lat,
output = "points"
) %>%
bind_cols() %>%
bind_cols(grid) %>%
mutate(year = yr,
region = reg,
spp = taxa) %>%
select(index, state, year, region, spp, lon=x, lat=y, value = z)
}
return(fit_tin)
}
# Test me
# Needs variables in Control panel
# Test no data: "Alosa aestivalis", reg = "Northeast US Fall", yr = 1974
# tis(input_data = ocean_data, grid = grid_eez_df, taxa = "Illex illecebrosus", reg = "Northeast US Fall", yr = 1973)
# lapply(years,tis,input_data = ocean_data, grid = grid_eez_df, taxa = "Illex illecebrosus", reg = regions[2])
Load required data. Note that some of it has been previously created
head(spp)
[1] "Alosa aestivalis" "Alosa pseudoharengus"
[3] "Alosa sapidissima" "Amblyraja radiata"
[5] "Ammodytes dubius" "Anarhichas lupus"
# single species run
# run_tis(input_data = ocean_data,
# grid = grid_eez_df,
# years = years,
# reg = regions,
# taxa = "Illex illecebrosus"
# )
# Run them all in parallel
parallel::mclapply(spp,
run_tis,
input_data = ocean_data,
grid = grid_eez_df,
years = years,
reg = regions
)
This graph shows the proportion of each State over time
This graph shows the proportion of each State smoothed over a 10 years running mean. It helps seeing the trends better (I Think…)